Introduction

This analysis has been performed to accompany a tomato growth trial in the Hydroponic Lab at the KSU Field Station (coordinates?). The trial was conducted with the purpose of determining the effects of differing soil microbial inoculation timings on tomato plant photosynthesis and crop yield and quality. Photosynthesis metrics were taken with a Li-COR Li-600 and two PhotosynQ MultispeQ v2.0s, and include stomatal conductance (gsw), measured in (mol m-2 s-1), which is a measure of how well the leaf is able to conduct gas through its stomates, and the efficiency of photosystem II (PhiPS2), where photosystem II is a piece of machinery in the chloroplast essential to photosynthesis and PhiPS2 is the quantum yield in light calculated from fluorescence.

Sample size (n=48) with 4 groups, 12 replicates per group. The tomato plants were grown in a dutch bucket hydroponic system from April to October 2024. Fluorometric measurements were taken twice weekly, and fruit were harvested upong fully ripening and taken back to the lab on the KSU campus for analysis.

Fruit analysis included weighing the tomato for its mass in grams and visually assessing it for blossom end rot (BER), a nutritional disorder caused by calcium deficiency that renders the fruit unmarketable, as well as visual checks for fungus and cracking. Fungus commonly appears on parts of the tomato with blossom end rot or cracking. From there, fruit were analyzed with a penetrometer for ripeness, (1-kg) where kg is mapped from 0:1, and then with Fisher Brix Refractometers for sugar concentration in percent sugar content of the tomato’s juice.

Discarded fruit were then composted.


1. Setup

1.1 Load packages

DATA WRANGLING - tidyverse

STATISTICS - lme4 - car - MASS - rstatix - pwr

GRAPHING - showtext - scico - ggpubr

OTHER - devtools - multiFitR - multiFitRgg - ztils

1.2 Graphical Setup

Load colors. Set the palette to one of the scientific color map palettes (scico) and create variants that contain the specific values for 2, 4, and 10 colors. This helps when working with factors, as it can be a pain to apply gradients otherwise. I also make the vector of colors =n+1 because the final value in scico palettes are often white, rendering them hard to see against a white plot background. (and we use white plot backgrounds because colored backgrounds are difficult to parse and expensive to print)

Load fonts. Google’s Open Sans for titles and Montserrat for body text.

Set shape selection from ggplot’s shape library.

a_palette <- "bilbao"
# super lazy way to make n+1 palettes
two_colors = scico(3, palette=a_palette)
four_colors = scico(5, palette=a_palette)
five_colors = scico(6, palette=a_palette)
true_two_col = scico(2, palette=a_palette)
ten_col = scico(10, palette=a_palette)
twelve_col = scico(13, palette=a_palette)

# add fonts
font_add_google("Open Sans", family = "open")
font_add_google("Montserrat", family = "mont")
# necessary to initialize imported fonts
showtext_auto()

# custom shapes
four_shapes = c(15,16,17,23)

1.3 Experimental Constants

## germination date
germdate <- "2024-05-01"

## treatment order
treatment_order <- c("Control",
                     "Transplantation",
                     "Germination", 
                     "Germ+Trans")

## plant factor order
plant_order <- c("1 1", "1 2", "1 3", "1 4", "1 5", "1 6",
                 "1 7", "1 8", "1 9", "1 10", "1 11", "1 12",
                 "2 1", "2 2", "2 3", "2 4", "2 5", "2 6",
                 "2 7", "2 8", "2 9", "2 10", "2 11", "2 12",
                 "3 1", "3 2", "3 3", "3 4", "3 5", "3 6",
                 "3 7", "3 8", "3 9", "3 10", "3 11", "3 12",
                 "4 1", "4 2", "4 3", "4 4", "4 5", "4 6",
                 "4 7", "4 8", "4 9", "4 10", "4 11", "4 12")

2. Load Data

2.1 Photosynthesis data

2.1.1 Li-600 data

Load Li-600 data from file and convert strings to factors. Then change the treatments from numbers to the names of the groups. From here, filter the data for only those readings that have a leak percent of less than 10. QUESTION: What leak percent is acceptable to filter out? Rename the Date column to Date_ref to use it as a check that we formatted dates correctly. Format date to POSIXct in month-day-year format. Create a new plant factor that corresponds to each unique plant, then re-level the Treatment for consistency across analyses. Then we make a new column that is the log of gsw.

From here, outliers in the main target variables (stomatal conductance [gsw] and photosystem II efficiency [PhiPS2]) were identified.

This is one of those things I might not bother with, as it is almost certainly not necessary to remove outliers.

Then, new dataframes were created out of the base Li-600 dataset to look at the mean and sd of gsw and PhiPS2, grouped by 1. treatment and date and 2. plant factor.

Another dataframe was created looking at mass and fruit count sums by treatment and plant.

As a final note, this data has something wrong with it, as there are ~100 gsw values <0, which should be impossible. We’ve contacted Li-COR for clarification and hopefully there’s just a correction we can apply to the data. If not, then Li-COR has some explaining to do.

## Load and clean Li-600 data
Li_data_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_LI600.csv"
Li_data <- read.csv(Li_data_file, stringsAsFactors = T) %>%
  mutate(Treatment = case_when(
          Row==1~"Control",
          Row==2~"Transplantation",
          Row==3~"Germination",
          Row==4~"Germ+Trans",
          TRUE~NA)) %>%
  filter(leak_pct<10, gsw>0) %>%
  rename(Date_ref = Date) %>%
  mutate(Date = parse_date_time(Date_ref, orders = "mdy"),
         Time = parse_date_time(Time, orders = "T"),
         DaysFromGermination = as.numeric(round(difftime(Date, germdate, units = c("days")), 0)),
         plant_fac = as.factor(paste(Row, Pot)),
         Treatment = factor(Treatment, levels = treatment_order),
         logitPS2 = logit(PhiPS2, FALSE),
         plant_fac = factor(plant_fac, levels = plant_order)
) %>%
  group_by(DaysFromGermination) %>%
  mutate(MinutesFromStart = round(difftime(Time, min(Time), units = "mins"), 4)) %>%
  ungroup() %>%
  mutate(
    Time = format(Time, "%H:%M:%S")
  )

# New dataframes of Li_data mean and sd for stomatal conductance (gsw) and Photosystem II efficiency (PhiPS2), grouped by treatment and plant
Li_data_stats_byplant <- Li_data %>%
  group_by(Treatment, plant_fac) %>%
  summarise_at(vars(gsw, logitPS2), list(mean=mean, sd=sd))

2.1.2 MultispeQ data

mq_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_Multispeq.csv"
m_data <- read.csv(mq_file) %>%
  mutate(Treatment = case_when(
      Row=="A"~"Control",
      Row=="B"~"Transplantation",
      Row=="C"~"Germination",
      Row=="D"~"Germ+Trans",
      TRUE~NA),
    Row_num = case_when(
      Row=="A"~ 1,
      Row=="B"~ 2,
      Row=="C"~ 3,
      Row=="D"~ 4,
      TRUE~NA),
    plant_fac = as.factor(paste(Row_num, Pot)),
    Date = as.Date(time, "%m/%d/%Y"),
    Time = parse_date_time(time, "%m/%d/%Y %H:%M %p"),
    DaysFromGermination = as.numeric(round(difftime(Date, germdate, units = c("days")), 0)),
    Treatment = factor(Treatment, levels = treatment_order),
     plant_fac = factor(plant_fac, levels = plant_order),
    logitPS2 = logit(Phi2, FALSE),
    logitFvPFmP = logit(FvP_over_FmP, FALSE)
    ) %>%
  group_by(DaysFromGermination) %>%
  mutate(MinutesFromStart = round(difftime(Time, min(Time), units = "mins"), 4)) %>%
  ungroup() %>%
  mutate(
    Time = format(Time, "%H:%M:%S")
  )

m_data_stats_byplant <- m_data %>%
  group_by(Treatment, plant_fac) %>%
  summarise_at(vars(logitFvPFmP, logitPS2), list(mean=mean, sd=sd))

2.2 Fruit Data

Load fruit data from set file name and location, converting strings to factors and renaming Treatment’s group names. Then setup “fruit” as a way to sum the total fruit count, create references for our dates then parse them to POSIXct in month-day-year format, then grab just the days and calculate the difference in days from harvest to analysis. (NOTE: this is susceptible to month rollover errors, so be sure not to make any of those mistakes) Make a plant factor as a unique ID for each plant, relevel the treatment for consistency across analyses, create mass bins (n=10) to look at BER occurrence across mass groups, and convert BER, fungus, and cracking to factors. Then, we map the penetrometer data from 0:1 and subtract it from 1.0, effectively giving us the ripeness. (this is because higher penetrometer values correspond to lower ripeness)

Then create a copy of the data that doesn’t have BER and another copy of just those with BER.

Make a pivot table of the means and sd of the tomato mass (g) and sugar (%), grouped by treatment.

Make another pivot table to get the sum of fruit count, mass, and BER, fungus, and cracking occurrence by treatment group. Then calculate the probability of BER, fungus, and cracking for each group. Do that again grouping the fruit into mass bins.

## Load and clean TIP24 fruit data
Fl_data_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_Fruit.csv"
Fl_data <- read.csv(Fl_data_file, stringsAsFactors = T) %>%
  filter(mass >0) %>%
  mutate(Treatment = case_when(
    row==1~"Control",
    row==2~"Transplantation",
    row==3~"Germination",
    row==4~"Germ+Trans",
    TRUE~NA),
         fruit = 1,
         date_analysis_ref = date_analysis,
         date_harvest_ref = date_harvest,
         date_analysis = parse_date_time(date_analysis_ref, orders = "mdy"),
         date_harvest = parse_date_time(date_harvest_ref, orders="mdy"),
         DaysFromHarvestToAnalysis = as.numeric(round(difftime(date_analysis,
                                                               date_harvest, units = c("days")), 0)),
         DaysFromGermination = as.numeric(round(difftime(date_analysis,
                                                         germdate, units = c("days")), 0)),
         plant_fac = as.factor(paste(row, plant)),
         Treatment = factor(Treatment, levels = treatment_order),
         plant_fac = factor(plant_fac, levels = plant_order),
         mass_bin = cut(mass, breaks=10),
         BER_fac = as.factor(BER),
         fungus_fac = as.factor(fungus),
         cracking_fac = as.factor(cracking),
         ripeness = abs(1 - round(penetrometer/max(na.omit(penetrometer)), 2))
         )

## Filter to only rows where there is no BER and sugar is greater than 0
Fl_data_no_BER <- Fl_data %>%
  filter(BER==0 & sugar_avg > 0)

## New dataframes, filtering to only those where BER is equal to 1
Fl_data_BER <- Fl_data %>%
  filter(BER==1)

## Make summary dataframe with mean and sd of mass and sugar by plant.
Fl_means_byplant <- Fl_data_no_BER %>%
  group_by(Treatment, plant_fac) %>%
  summarise_at(vars(mass, sugar_avg), list(mean=mean, sd=sd))

## Sum the fruit, BER, fungus, cracking, and mass by treatment group, then get probs of each factor.
Fl_data_summary <- Fl_data %>%
  group_by(Treatment) %>%
  summarise_at(vars(fruit, mass, BER, fungus, cracking, mass),
               list(sum=sum)) %>%
  mutate(pBER = round(BER_sum/fruit_sum, 4),
         pfungus = round(fungus_sum/fruit_sum, 4),
         pcracking = round(cracking_sum/fruit_sum, 4)
         )

## Sum the fruit, BER, fungus, and cracking by mass bin, then get probs of each factor.
Fl_data_mb <- Fl_data %>%
  group_by(Treatment, mass_bin) %>%
    summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
  mutate(pBER = round(BER_sum/fruit_sum, 4),
         pfungus = round(fungus_sum/fruit_sum, 4),
         pcracking = round(cracking_sum/fruit_sum, 4)
         )

## Sum the fruit, BER, fungus, and cracking by plant, then get probs of each factor.
Fl_sum_byplant <- Fl_data %>%
  group_by(Treatment, plant) %>%
    summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
  mutate(pBER = round(BER_sum/fruit_sum, 4),
         pfungus = round(fungus_sum/fruit_sum, 4),
         pcracking = round(cracking_sum/fruit_sum, 4)
         )

Fl_sum_means <- Fl_sum_byplant %>%
  group_by(Treatment) %>%
  summarise_at(vars(fruit_sum, mass_sum, pBER), list(mean=mean, sd=sd))

Fl_means <- Fl_data %>%
  group_by(Treatment) %>%
  summarise_at(vars(mass), list(mean=mean, sd=sd))

3. Analyses

3.1 Photosynthesis Analysis

3.1.1 Stomatal Conductance (GSW) exploratory

# PDF, CDF, and KS test
multiPDF_Z(Li_data$gsw, 100, "all", a_palette, "gsw")

multiCDF_Z(Li_data$gsw, 100, "all", a_palette, "gsw")

multiKS_cont(Li_data$gsw, "all")
##   Distribution Distance PValue
## 1       Normal    0.165  0.000
## 2    Lognormal    0.089  0.000
## 3        Gamma    0.058  0.039
## 4  Exponential    0.139  0.000
# check for heteroscedasticity
## in gsw as a function of treatment
leveneTest(Li_data$gsw~Li_data$Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   3  2.8428 0.03717 *
##       578                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## in gsw as a function of time from start
leveneTest(Li_data$gsw~Li_data$MinutesFromStart)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group 342  0.9056 0.7996
##       239
## gsw by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = gsw, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.15, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("Stomatal Conductance Across Inoculation Treatments in Tomato", 40)
  ) +
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

# Mean GSW per plant by Treatment boxplot 
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = gsw_mean, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.15, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("Mean Stomatal Conductance Per Plant
                   Across Inoculation Treatments in Tomato", 50)
  ) +
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

# GSW by date across relative humidity
ggplot(data = Li_data, aes(x=DaysFromGermination, y = gsw, 
                           color = rh_s)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_scico(begin=0.9, end=0, palette=a_palette)+
#  scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
  theme_bw()+
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  xlab("Days From Germination")+
  labs(
    title=str_wrap("Stomatal Conductance By Days From Germination Across Relative Humidity in Tomato", 50)
  )+
    guides(color=guide_colorbar(title=str_wrap("Relative Humidity", 10)))+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.position="right",
    legend.title.position = "top",
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=two_colors[2], fill=NA,
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

# GSW by time across relative humidity
ggplot(data = Li_data, aes(x=Time, y = gsw, 
                           color = Date_ref,
                           shape = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_shape_manual(values=four_shapes)+
  scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
  theme_bw()+
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  xlab("Time")+
  labs(
    title=str_wrap("Stomatal Conductance By Time Across Date in Tomato", 50)
  )+
    guides(color=guide_legend(title="Date"))+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.position="right",
    legend.title.position = "top",
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=two_colors[2], fill=NA,
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

# gsw by relative humidity across treatment groups
ggplot(data = Li_data, aes(x=rh_s, y = gsw, 
                           shape=Treatment, fill=Treatment, color = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_fill_manual(values=four_colors)+
  #ylim(0,1)+
  theme_bw()+
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  xlab("Relative Humidity")+
  labs(
    title=str_wrap("Stomatal Conductance By Relative
                   Humidity Across Inoculation Treatment in Tomato", 50)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.position="inside",
    legend.title.position = "top",
    legend.position.inside=c(0.15,0.75),
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

# gsw by Minutes from Start
ggplot(data = Li_data, aes(x=MinutesFromStart, y = gsw, color = rh_s)) +
  geom_point()+
  scale_color_scico(begin=0.9, end=0, palette=a_palette)+
  facet_wrap(~Treatment)+
  theme_bw()+
  ylab("Stomatal Conductance (mol m-2 s-1)")+
  xlab("Minutes From Start")+
  labs(
    title=str_wrap("Stomatal Conductance by Minutes From Start Across Inoculation Treatment in Tomato", 50)
  )+
      guides(color=guide_colorbar(title=str_wrap("Relative Humidity", 8)))+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.

3.1.3 Stomatal Conductance (GSW) Analysis

Do a generalized linear mixed model with stomatal conductance (gsw) as response, Treatment as predictor, and days from germination as a fixed effect that should incorporate fluctuations in environmental variables, such as relative humidity, with a gamma family, log link function because gsw best fits a gamma distribution.

#GLMM with Treatment as explanatory and days from germination as fixed effect
gsw_glmm_dfg_summary <- summary(gsw_glm_dfg <- (glm(gsw ~ Treatment + DaysFromGermination, 
                  data=Li_data, family=Gamma(link="log"))))
print(gsw_glmm_dfg_summary)
## 
## Call:
## glm(formula = gsw ~ Treatment + DaysFromGermination, family = Gamma(link = "log"), 
##     data = Li_data)
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -0.545480   0.120872  -4.513 7.75e-06 ***
## TreatmentTransplantation -0.215724   0.099400  -2.170  0.03039 *  
## TreatmentGermination     -0.285432   0.100647  -2.836  0.00473 ** 
## TreatmentGerm+Trans      -0.541036   0.102431  -5.282 1.81e-07 ***
## DaysFromGermination       0.003901   0.002176   1.792  0.07361 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.7609672)
## 
##     Null deviance: 418.91  on 581  degrees of freedom
## Residual deviance: 395.12  on 577  degrees of freedom
## AIC: 396.87
## 
## Number of Fisher Scoring iterations: 6
#GLMM with Treatment as explanatory and relative humidity as fixed effect
gsw_glmm_rhs_summary <- summary(gsw_glm_rhs <- (glm(gsw ~ Treatment + rh_s, 
                  data=Li_data, family=Gamma(link="log"))))
print(gsw_glmm_rhs_summary)
## 
## Call:
## glm(formula = gsw ~ Treatment + rh_s, family = Gamma(link = "log"), 
##     data = Li_data)
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -5.903211   0.199662 -29.566  < 2e-16 ***
## TreatmentTransplantation -0.248352   0.057316  -4.333 1.73e-05 ***
## TreatmentGermination     -0.269800   0.058169  -4.638 4.35e-06 ***
## TreatmentGerm+Trans      -0.269977   0.059667  -4.525 7.35e-06 ***
## rh_s                      0.071222   0.002588  27.519  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.2532233)
## 
##     Null deviance: 418.91  on 581  degrees of freedom
## Residual deviance: 182.51  on 577  degrees of freedom
## AIC: -87.276
## 
## Number of Fisher Scoring iterations: 9
# pseudo R2 of glmms
gsw_glmm_dfg_pseudoR2 <- (gsw_glmm_dfg_summary$null.deviance -
                        gsw_glmm_dfg_summary$deviance)/gsw_glmm_dfg_summary$null.deviance
print(gsw_glmm_dfg_pseudoR2)
## [1] 0.05678831
gsw_glmm_rhs_pseudoR2 <- (gsw_glmm_rhs_summary$null.deviance -
                        gsw_glmm_rhs_summary$deviance)/gsw_glmm_rhs_summary$null.deviance
print(gsw_glmm_rhs_pseudoR2)
## [1] 0.5643147
# PREDICTIONS
  ## im going to do this in maybe the worst possible way, sorry

# GSW_RHS
## initialize prediction plot
p <- ggplot()+
  geom_point(data=Li_data, aes(x=rh_s, y=gsw, color=Treatment, shape=Treatment),
             alpha = 0.8)

## base GLM
gsw_rhs <- (glm(gsw ~ rh_s, 
                data=Li_data, 
                family=Gamma(link="log")))
## sequence
rh_s <- seq(min(Li_data$rh_s), max(Li_data$rh_s), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  gsw_rhs <- (glm(gsw ~ rh_s, 
                data=Li_data %>% filter(Treatment == x), #treatmnent == x
                family=Gamma(link="log")))
  ## calculate predictions
  prx <- as.data.frame(x=rh_s)
  pred <- predict(gsw_rhs, newdata=prx, se.fit=TRUE, type="link")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=rh_s, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=rh_s, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=rh_s, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=rh_s, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
  facet_wrap(~Treatment)+
  labs(title = "Real vs Predicted for gsw_rhs glm")+
  ylab("Stomatal Conductance (mol+1 m-2 s-1)")+
  xlab("Relative Humidity (%)")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="none",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

# GSW_DFG

## initialize prediction plot
p <- ggplot()+
  geom_point(data=Li_data, aes(x=DaysFromGermination, y=gsw, color=Treatment, shape=Treatment),
             alpha = 0.8)

## base GLM
gsw_dfg <- (glm(gsw ~ DaysFromGermination, 
                data=Li_data, 
                family=Gamma(link="log")))
## sequence
DaysFromGermination <- seq(min(Li_data$DaysFromGermination), max(Li_data$DaysFromGermination), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  glm_x <- (glm(gsw ~ DaysFromGermination, 
                data=Li_data %>% filter(Treatment == x), #treatmnent == x
                family=Gamma(link="log")))
  ## calculate predictions
  prx <- as.data.frame(x=DaysFromGermination)
  pred <- predict(glm_x, newdata=prx, se.fit=TRUE, type="link")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
  facet_wrap(~Treatment)+
  labs(title = "Real vs Predicted for gsw_dfg glm")+
  ylab("Stomatal Conductance (mol+1 m-2 s-1)")+
  xlab("Days From Germination")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="none",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

3.1.3 PhiPS2 exploratory

# Li-600
## PDF and CDF plots
multiPDF_Z(subset(Li_data$logitPS2, Li_data$logitPS2>0), 100, "all", a_palette, "Logit PhiPS2")

multiCDF_Z(subset(Li_data$logitPS2, Li_data$logitPS2>0), 100, "all", a_palette, "Logit PhiPS2")

## multi KS tests
multiKS_cont(subset(Li_data$logitPS2, Li_data$logitPS2>0), "all")
## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test
##   Distribution Distance PValue
## 1       Normal    0.310      0
## 2    Lognormal    0.157      0
## 3        Gamma    0.212      0
## 4  Exponential    0.449      0
## test for heteroscedasticity
leveneTest(Li_data$logitPS2~Li_data$Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   3   2.036 0.1077
##       578
leveneTest(Li_data$logitPS2~Li_data$MinutesFromStart)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group 342  0.7906 0.9766
##       239
# Li PhiPS2 by Days From Germination, colored by treatment
ggplot(data = Li_data, aes(x=DaysFromGermination, y = PhiPS2, color = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  ylim(.6,.8)+
  theme_bw()+
  ylab("PhiPS2")+
  xlab("Days From Germination")+
  labs(
    title=str_wrap("Li-600 Photosystem II Efficiency by Date Across Inoculation Treatment in Tomato", 50)
  )+
      guides(color=guide_legend(title=str_wrap("Treatment", 8)))+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

# PhiPS2 by Minutes from Start colored by treatment
ggplot(data = Li_data, aes(x=MinutesFromStart, y = PhiPS2, color = Treatment)) +
  geom_point()+
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  theme_bw()+
  ylab("PhiPS2")+
  xlab("Minutes From Start")+
  labs(
    title=str_wrap("Li-600 Photosystem II Efficiency by Minutes From Start Across Inoculation Treatment in Tomato", 50)
  )+
      guides(color=guide_legend(title=str_wrap("Treatment", 8)))+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.

## logitPS2 by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = logitPS2, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.05, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("Li-600 LogitPS2 Across Inoculation Treatments in Tomato", 40)
  ) +
  ylab("Logit PhiPS2")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

# Mean logitPS2 per plant by Treatment boxplot 
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = logitPS2_mean, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.15, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("Li-600 Mean Logit PhiPS2 Per Plant
                   Across Inoculation Treatments in Tomato", 60)
  ) +
  ylab("Logit PhiPS2")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

# MultispeQ
## PDF, CDF, and KS test
multiPDF_Z(subset(m_data$logitPS2, m_data$logitPS2>0),
              100, c('normal','lognormal','gamma'), a_palette, "Logit PhiPS2")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

multiCDF_Z(subset(m_data$logitPS2, m_data$logitPS2>0),
              100, c('normal','lognormal','gamma'), a_palette, "Logit PhiPS2")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

multiKS_cont(subset(m_data$logitPS2, m_data$logitPS2>0),
             c('normal','lognormal','gamma'))
## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
##   Distribution Distance PValue
## 1       Normal    0.048  0.070
## 2    Lognormal    0.052  0.043
## 3        Gamma    0.034  0.393
# check for heteroscedasticity
leveneTest(m_data$logitPS2~m_data$Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   3  3.6607 0.01224 *
##       718                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(data = m_data, aes(x=DaysFromGermination, y = logitPS2, 
                           shape=Treatment, fill=Treatment, color = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_fill_manual(values=four_colors)+
  theme_bw()+
  ylab("Logit PhiPS2")+
  xlab("Days From Germination")+
  labs(
    title=str_wrap("MultispeQ Logit Photosystem II Efficiency By Date Across
                   Inoculation Treatment in Tomato", 50)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

ggplot(data = m_data, aes(x=MinutesFromStart, y = logitPS2, 
                           shape=Treatment, fill=Treatment, color = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_fill_manual(values=four_colors)+
  theme_bw()+
  ylab("Logit PhiPS2")+
  xlab("Minutes From Start")+
  labs(
    title=str_wrap("MultispeQ Logit Photosystem II Efficiency By Minutes From Start Across
                   Inoculation Treatment in Tomato", 50)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.

## logitPS2 by Treatment boxplot
ggplot(data = m_data, aes(x= Treatment, y = logitPS2, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.15, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("MultispeQ LogitPS2 Across Inoculation Treatments in Tomato", 40)
  ) +
  ylab("Logit PhiPS2")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

# Mean logitPS2 per plant by Treatment boxplot 
ggplot(data = m_data_stats_byplant, aes(x= Treatment, y = logitPS2_mean, fill=Treatment, color=Treatment)) +
  geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
    geom_jitter( width=.15, height=0)+
#    stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
#  stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  labs(
    title=str_wrap("MultispeQ Mean Logit PhiPS2 Per Plant
                   Across Inoculation Treatments in Tomato", 50)
  ) +
  ylab("Mean Logit PhiPS2")+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## PhiPS2 by date, grouped by treatment
ggplot(data = Li_data, aes(x=as.factor(Date), y = PhiPS2, color = Treatment,
                           fill=Treatment, shape=Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_fill_manual(values=four_colors)+
  scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
  ylim(.6,.8)+
  facet_wrap(~Treatment)+
  theme_bw()+
  ylab("PhiPS2")+
  xlab("Date")+
  labs(
    title=str_wrap("PhiPS2 By Date Across Microbial Treatments in Tomato", 60)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="none",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    axis.text.x = element_text(size=12, family = "mont", angle=10),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

# PhiPS2 by leaf across temperature
ggplot(data = Li_data, aes(x=Date, y = PhiPS2, color = Tleaf)) +
  geom_jitter(width=1, size=2, alpha=1, aes(x=as.factor(Date)))+
  scale_color_scico(begin=0.9, end=0, palette=a_palette)+
  scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
  ylim(.6,.8)+
  theme_bw()+
  guides(color=guide_legend(title="Leaf Temp (C)"))+
  ylab("PhiPS2")+
  xlab("Date")+
  labs(
    title=str_wrap("PhiPS2 By Date Across Temperature in Tomato", 50)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.position="inside",
    legend.title.position = "top",
    legend.position.inside=c(0.75,0.2),
    legend.key.height = unit(.3, "cm"),
   legend.background = element_rect(color=two_colors[2], 
                                    fill=NA, linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

3.1.4 PhiPS2 Analysis

Because we logit transformed PhiPS2, we should be able to treat it as normal.

# linear mixed model with PhiPS2 as response, Treatment as predictor (explanatory), days from germination as fixed effect

mod_Phi2_dfg_LI600_s <- summary(mod_Phi2_dfg_LI600 <- (lm(
  logitPS2 ~ Treatment + DaysFromGermination,
  data = Li_data)))
print(mod_Phi2_dfg_LI600_s)
## 
## Call:
## lm(formula = logitPS2 ~ Treatment + DaysFromGermination, data = Li_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9216 -0.1892 -0.0868  0.0226  4.4045 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.477569   0.079358   6.018 3.14e-09 ***
## TreatmentTransplantation 0.158414   0.065261   2.427  0.01551 *  
## TreatmentGermination     0.097558   0.066080   1.476  0.14039    
## TreatmentGerm+Trans      0.203130   0.067251   3.020  0.00264 ** 
## DaysFromGermination      0.008719   0.001429   6.102 1.93e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5727 on 577 degrees of freedom
## Multiple R-squared:  0.07766,    Adjusted R-squared:  0.07126 
## F-statistic: 12.15 on 4 and 577 DF,  p-value: 1.741e-09
## R2 0.0713

mod_Phi2_qamb_LI600_s <- summary(mod_Phi2_dfg_LI600 <- (lm(
  logitPS2 ~ Treatment + Qamb,
  data = Li_data)))
print(mod_Phi2_qamb_LI600_s)
## 
## Call:
## lm(formula = logitPS2 ~ Treatment + Qamb, data = Li_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8591 -0.1740 -0.0763  0.0093  4.6414 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.0644152  0.0889739  11.963  < 2e-16 ***
## TreatmentTransplantation  0.2112193  0.0687409   3.073  0.00222 ** 
## TreatmentGermination      0.0784828  0.0688119   1.141  0.25453    
## TreatmentGerm+Trans       0.2283420  0.0695834   3.282  0.00109 ** 
## Qamb                     -0.0006152  0.0002478  -2.483  0.01333 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5878 on 577 degrees of freedom
## Multiple R-squared:  0.02852,    Adjusted R-squared:  0.02179 
## F-statistic: 4.235 on 4 and 577 DF,  p-value: 0.002186
## R2 0.0218

mod_Phi2_dfg_Mult_s <- summary(mod_Phi2_dfg_Mult <- (lm(
  logitPS2 ~ Treatment + DaysFromGermination,
  data = m_data)))
print(mod_Phi2_dfg_Mult_s)
## 
## Call:
## lm(formula = logitPS2 ~ Treatment + DaysFromGermination, data = m_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60212 -0.10424 -0.00623  0.09618  0.46918 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.4217061  0.0228180  18.481  < 2e-16 ***
## TreatmentTransplantation  0.0968281  0.0155648   6.221  8.4e-10 ***
## TreatmentGermination     -0.0394771  0.0154092  -2.562   0.0106 *  
## TreatmentGerm+Trans       0.0394109  0.0156349   2.521   0.0119 *  
## DaysFromGermination       0.0004075  0.0002951   1.381   0.1677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1488 on 717 degrees of freedom
## Multiple R-squared:  0.1048, Adjusted R-squared:  0.09982 
## F-statistic: 20.99 on 4 and 717 DF,  p-value: 2.223e-16
## R2 0.0998

mod_Phi2_qamb_Mult_s <- summary(mod_Phi2_qamb_Mult <- (lm(
  logitPS2 ~ Treatment + Light.Intensity..PAR.,
  data = m_data)))
print(mod_Phi2_qamb_Mult_s)
## 
## Call:
## lm(formula = logitPS2 ~ Treatment + Light.Intensity..PAR., data = m_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40277 -0.06556 -0.00085  0.07001  0.32258 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               7.129e-01  1.344e-02  53.039  < 2e-16 ***
## TreatmentTransplantation  3.563e-02  1.179e-02   3.021  0.00261 ** 
## TreatmentGermination     -8.851e-02  1.159e-02  -7.639    7e-14 ***
## TreatmentGerm+Trans      -2.665e-02  1.189e-02  -2.241  0.02530 *  
## Light.Intensity..PAR.    -9.182e-04  3.767e-05 -24.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1102 on 717 degrees of freedom
## Multiple R-squared:  0.5092, Adjusted R-squared:  0.5065 
## F-statistic:   186 on 4 and 717 DF,  p-value: < 2.2e-16
## R2 0.5065


# run and plot predictions

# Li-600 logitPS2 ~ Qamb
## initialize prediction plot
p <- ggplot()+
  geom_point(data=Li_data, aes(x=Qamb, y=logitPS2, color=Treatment, shape=Treatment),
             alpha = 0.8)

## base LM
lm_x <- (lm(logitPS2 ~ Qamb, data=Li_data))
## sequence
Qamb <- seq(min(Li_data$Qamb), max(Li_data$Qamb), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  lm_x <- (lm(logitPS2 ~ Qamb, 
                data=Li_data %>% filter(Treatment == x)))
  ## calculate predictions
  prx <- as.data.frame(x=Qamb)
  pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=Qamb, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=Qamb, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=Qamb, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=Qamb, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for Li-600 logitPS2 ~ Treatment + Qamb")+
  ylab("logit PhiPS2")+
  xlab("Ambient Light (lumens)")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

# Li-600 logitPS2 ~ DaysFromGermination
## initialize prediction plot
p <- ggplot()+
  geom_point(data=Li_data, aes(x=DaysFromGermination, y=logitPS2, color=Treatment, shape=Treatment),
             alpha = 0.8)

## base LM
lm_x <- (lm(logitPS2 ~ DaysFromGermination, 
                data=Li_data))
## sequence
DaysFromGermination <- seq(min(Li_data$DaysFromGermination), max(Li_data$DaysFromGermination), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  lm_x <- (lm(logitPS2 ~ DaysFromGermination, 
                data=Li_data %>% filter(Treatment == x)))
  ## calculate predictions
  prx <- as.data.frame(x=DaysFromGermination)
  pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for Li-600 logitPS2 ~ Treatment + DaysFromGermination")+
  ylab("logit PhiPS2")+
  xlab("Days From Germination")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

# Multispeq logitPS2 ~ DaysFromGermination
## initialize prediction plot
p <- ggplot()+
  geom_point(data=m_data, aes(x=DaysFromGermination, y=logitPS2, color=Treatment, shape=Treatment),
             alpha = 0.8)
## base GLM
lm_x <- (lm(logitPS2 ~ DaysFromGermination, 
                data=m_data))
## sequence
DaysFromGermination <- seq(min(m_data$DaysFromGermination), max(m_data$DaysFromGermination), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  lm_x <- (lm(logitPS2 ~ DaysFromGermination, 
                data=m_data %>% filter(Treatment == x))) #treatmnent == x
  ## calculate predictions
  prx <- as.data.frame(x=DaysFromGermination)
  pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for MultispeQ logitPS2 ~ Treatment + DaysFromGermination")+
  ylab("logit PhiPS2")+
  xlab("Days From Germination")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

# Multispeq logitPS2 ~ Light.Intensity..PAR.
## initialize prediction plot
p <- ggplot()+
  geom_point(data=m_data, aes(x=Light.Intensity..PAR., y=logitPS2, color=Treatment, shape=Treatment),
             alpha = 0.8)
## base GLM
lm_x <- (lm(logitPS2 ~ Light.Intensity..PAR., 
                data=m_data))
## sequence
Light.Intensity..PAR. <- seq(min(m_data$Light.Intensity..PAR.), max(m_data$Light.Intensity..PAR.), length=200)

# loop over treatments
for (x in treatment_order) {
 ## treatment specific GLM
  lm_x <- (lm(logitPS2 ~ Light.Intensity..PAR., 
                data=m_data %>% filter(Treatment == x))) #treatmnent == x
  ## calculate predictions
  prx <- as.data.frame(x=Light.Intensity..PAR.)
  pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
  prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
  prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
  prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
  prx$Treatment <- x
  prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
  ## add predicted points
  p <- p + 
  geom_line(data=prx, aes(x=Light.Intensity..PAR., y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_line(data=prx, aes(x=Light.Intensity..PAR., y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
  geom_line(data=prx, aes(x=Light.Intensity..PAR., y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
  geom_ribbon(data=prx, aes(x=Light.Intensity..PAR., ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
  scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
  scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for MultispeQ logitPS2 ~ Treatment + Light")+
  ylab("logit PhiPS2")+
  xlab("Light Intensity (PAR)")+
  theme_minimal()+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position="right",
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )

print(p)

3.1.5 FvP/FmP exploratory

multiPDF_Z(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
              100, c('normal','lognormal','gamma'), a_palette, "FvP/FmP")

multiCDF_Z(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
              100, c('normal','lognormal','gamma'), a_palette, "FvP/FmP")

multiKS_cont(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
             c('normal','lognormal','gamma'))
## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
##   Distribution Distance PValue
## 1       Normal    0.068  0.002
## 2    Lognormal    0.057  0.018
## 3        Gamma    0.058  0.015
leveneTest(m_data$logitFvPFmP~m_data$Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   3   1.595 0.1892
##       718
ggplot(data = m_data, aes(x=MinutesFromStart, y = logitFvPFmP, 
                           shape=Treatment, fill=Treatment, color = Treatment)) +
  geom_jitter(width=1, size=2, alpha=1)+
  scale_color_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_fill_manual(values=four_colors)+
  facet_wrap(~Treatment)+
  theme_bw()+
  ylab("Logit FvP/FmP")+
  xlab("Minutes From Start")+
  labs(
    title=str_wrap("Logit FvP/FmP By Minutes From Start Across
                   Inoculation Treatment in Tomato", 50)
  )+
  theme(
    text = element_text(size=24, family="mont"),
    legend.position = "none",
    legend.title = element_text(size=24, family="mont", face="bold"),
    legend.text = element_text(size=20, family="mont"),
    legend.key.height = unit(.3, "cm"),
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold", lineheight = .5)
  )
## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.

3.1.6 FvP/FmP analysis

Because FvP/FmP is a unitless ratio from 0:1, we can logit transform it and then treat it as if it were normal.

# Linear model of logit transformed FvPFmP as a function of Treatment with days from germination as a fixed effect
mod_FvP_dfg_s <- summary(mod_FvP_dfg <- (lm(
  logitFvPFmP ~ Treatment + DaysFromGermination,
  data = m_data)))
print(mod_FvP_dfg_s)
## 
## Call:
## lm(formula = logitFvPFmP ~ Treatment + DaysFromGermination, data = m_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.271209 -0.065279 -0.004234  0.071809  0.235475 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.1952892  0.0141777  84.308  < 2e-16 ***
## TreatmentTransplantation  0.0178982  0.0096710   1.851 0.064623 .  
## TreatmentGermination     -0.0371504  0.0095743  -3.880 0.000114 ***
## TreatmentGerm+Trans      -0.0226343  0.0097146  -2.330 0.020087 *  
## DaysFromGermination      -0.0007802  0.0001833  -4.256 2.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09244 on 717 degrees of freedom
## Multiple R-squared:  0.07253,    Adjusted R-squared:  0.06736 
## F-statistic: 14.02 on 4 and 717 DF,  p-value: 5.106e-11
## R2 0.06736

# Linear model of logit transformed FvPFmP as a function of Treatment with minutes from start as a fixed effect
mod_FvP_min_s <- summary(mod_FvP_min <- (lm(
  logitFvPFmP ~ Treatment + MinutesFromStart,
  data = m_data)))
print(mod_FvP_min_s)
## 
## Call:
## lm(formula = logitFvPFmP ~ Treatment + MinutesFromStart, data = m_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.281106 -0.064749 -0.008379  0.069584  0.231633 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.1716386  0.0081069 144.523  < 2e-16 ***
## TreatmentTransplantation  0.0645615  0.0120519   5.357 1.14e-07 ***
## TreatmentGermination     -0.0356205  0.0094415  -3.773 0.000175 ***
## TreatmentGerm+Trans       0.0267451  0.0124211   2.153 0.031636 *  
## MinutesFromStart         -0.0020292  0.0003227  -6.288 5.59e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09112 on 717 degrees of freedom
## Multiple R-squared:  0.0988, Adjusted R-squared:  0.09377 
## F-statistic: 19.65 on 4 and 717 DF,  p-value: 2.316e-15
## R2 0.09377

# Linear model of logit transformed FvPFmP as a function of Treatment with Light Intensity as a fixed effect
mod_FvP_l_s <- summary(mod_FvP_l<- (lm(
  logitFvPFmP ~ Treatment + Light.Intensity..PAR.,
  data = m_data)))
print(mod_FvP_l_s)
## 
## Call:
## lm(formula = logitFvPFmP ~ Treatment + Light.Intensity..PAR., 
##     data = m_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.250251 -0.070537 -0.007842  0.065663  0.264618 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.182e+00  1.127e-02 104.922  < 2e-16 ***
## TreatmentTransplantation  8.900e-03  9.886e-03   0.900 0.368287    
## TreatmentGermination     -4.478e-02  9.712e-03  -4.611 4.75e-06 ***
## TreatmentGerm+Trans      -3.309e-02  9.968e-03  -3.320 0.000947 ***
## Light.Intensity..PAR.    -1.398e-04  3.158e-05  -4.428 1.10e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09234 on 717 degrees of freedom
## Multiple R-squared:  0.07442,    Adjusted R-squared:  0.06926 
## F-statistic: 14.41 on 4 and 717 DF,  p-value: 2.521e-11
## R2 0.06926

# Linear model of logit transformed FvPFmP as a function of Treatment with Light Intensity  AND minutes from start as fixed effects
mod_FvP_lm_s <- summary(mod_FvP_lm<- (lm(
  logitFvPFmP ~ Treatment + Light.Intensity..PAR. + MinutesFromStart,
  data = m_data)))
print(mod_FvP_lm_s)
## 
## Call:
## lm(formula = logitFvPFmP ~ Treatment + Light.Intensity..PAR. + 
##     MinutesFromStart, data = m_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28136 -0.06668 -0.00711  0.06862  0.26022 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.218e+00  1.214e-02 100.298  < 2e-16 ***
## TreatmentTransplantation  5.698e-02  1.195e-02   4.770 2.23e-06 ***
## TreatmentGermination     -4.383e-02  9.425e-03  -4.650 3.95e-06 ***
## TreatmentGerm+Trans       1.848e-02  1.232e-02   1.500    0.134    
## Light.Intensity..PAR.    -1.553e-04  3.072e-05  -5.055 5.47e-07 ***
## MinutesFromStart         -2.149e-03  3.182e-04  -6.754 2.98e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0896 on 716 degrees of freedom
## Multiple R-squared:  0.1299, Adjusted R-squared:  0.1238 
## F-statistic: 21.37 on 5 and 716 DF,  p-value: < 2.2e-16
## R2 0.06926

3.2 Fruit

3.2.1 Mass exploratory

# PDF, CDF, and KS test
multiPDF_Z(Fl_data$mass, 100, "all", a_palette, "Mass")

multiCDF_Z(Fl_data$mass, 100, "all", a_palette, "Mass")

multiKS_cont(Fl_data$mass, "all")
## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test
##   Distribution Distance PValue
## 1       Normal    0.205  0.000
## 2    Lognormal    0.043  0.001
## 3        Gamma    0.104  0.000
## 4  Exponential    0.099  0.000
# Tests for homoscedasticity
leveneTest(Fl_data$mass~Fl_data$Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    3   5.993 0.0004599 ***
##       2080                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## plot of mass by treatment with box plots and jittered points
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=mass, color=Treatment, fill=Treatment))+
  geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  ylab("Tomato Mass (g)")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Tomato Mass by Treatment"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## plot of mass mean per plant by treatment with box plots and jittered points
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=mass_mean, color=Treatment, fill=Treatment))+
  geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  ylab("Mean Tomato Mass (g)")+
  xlab("Treatment")+
  theme_bw()+
  #facet_wrap(~Treatment)+
  labs(
    title=str_wrap("Mean Tomato Mass Per Plant across Treatments", 50)
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

3.2.2 Mass analysis

According to the PDFs, CDFs, and Kolmogorov-Smirnov tests, mass follows a Gamma or Lognormal distribution. I am unsure of whether to analyze it with each fruit as an observation or with each plant as an observation, where I’m looking at the mean tomato mass per plant instead of the individual tomato masses.

# INDIVIDUAL GAMMA glm
mT_glm_g_s <- summary(mT_glm_g <- glm(mass~Treatment, data=Fl_data, family=Gamma(link="log")))
print(mT_glm_g_s)
## 
## Call:
## glm(formula = mass ~ Treatment, family = Gamma(link = "log"), 
##     data = Fl_data)
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3.98752    0.04525  88.117  < 2e-16 ***
## TreatmentTransplantation  0.25348    0.06749   3.756 0.000178 ***
## TreatmentGermination      0.06497    0.06692   0.971 0.331704    
## TreatmentGerm+Trans      -0.04824    0.06490  -0.743 0.457365    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 1.183615)
## 
##     Null deviance: 1906.7  on 2083  degrees of freedom
## Residual deviance: 1879.7  on 2080  degrees of freedom
## AIC: 21005
## 
## Number of Fisher Scoring iterations: 5
mT_glm_g_pseudoR2 <- (mT_glm_g_s$null.deviance - mT_glm_g_s$deviance)/mT_glm_g_s$null.deviance
print(mT_glm_g_pseudoR2)
## [1] 0.0141688
## super high AIC and terrible PR2, so probably not the way to do it.

# GROUPED GAMMA glm
mTg_glm_g_s <- summary(mTg_glm <- glm(mass_mean~Treatment, data=Fl_means_byplant, family=Gamma(link="log")))
print(mTg_glm_g_s)
## 
## Call:
## glm(formula = mass_mean ~ Treatment, family = Gamma(link = "log"), 
##     data = Fl_means_byplant)
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               4.81769    0.07606  63.343   <2e-16 ***
## TreatmentTransplantation  0.07786    0.10756   0.724    0.473    
## TreatmentGermination      0.07615    0.10756   0.708    0.483    
## TreatmentGerm+Trans      -0.15292    0.10756  -1.422    0.162    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06941743)
## 
##     Null deviance: 3.3657  on 47  degrees of freedom
## Residual deviance: 2.9538  on 44  degrees of freedom
## AIC: 472.45
## 
## Number of Fisher Scoring iterations: 4
mTg_glm_g_pseudoR2 <- (mTg_glm_g_s$null.deviance - mTg_glm_g_s$deviance)/mTg_glm_g_s$null.deviance
print(mTg_glm_g_pseudoR2)
## [1] 0.1223773
# GROUPED Lognormal glm
mTg_glm_ln_s <- summary(mTg_glm_ln <- glm(mass_mean~Treatment, data=Fl_means_byplant, family=gaussian(link="log")))
print(mTg_glm_ln_s)
## 
## Call:
## glm(formula = mass_mean ~ Treatment, family = gaussian(link = "log"), 
##     data = Fl_means_byplant)
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               4.81769    0.07550  63.808   <2e-16 ***
## TreatmentTransplantation  0.07786    0.10286   0.757    0.453    
## TreatmentGermination      0.07615    0.10294   0.740    0.463    
## TreatmentGerm+Trans      -0.15292    0.11593  -1.319    0.194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1046.424)
## 
##     Null deviance: 52069  on 47  degrees of freedom
## Residual deviance: 46041  on 44  degrees of freedom
## AIC: 475.79
## 
## Number of Fisher Scoring iterations: 4
mTg_glm_ln_pseudoR2 <- (mTg_glm_ln_s$null.deviance - mTg_glm_ln_s$deviance)/mTg_glm_ln_s$null.deviance
print(mTg_glm_ln_pseudoR2)
## [1] 0.1157707
# Predictions
## How do I predict these?

# Post Hoc Tests

3.2.3 Sugar exploratory

# PDF, CDF, and KS test
multiPDF_Z(Fl_data_no_BER$sugar_avg, 100, "all", a_palette, "Sugar Concentration")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

multiCDF_Z(Fl_data_no_BER$sugar_avg, 100, "all", a_palette, "Sugar Concentration")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

multiKS_cont(Fl_data_no_BER$sugar_avg, "all")
## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test
##   Distribution Distance PValue
## 1       Normal    0.064  0.012
## 2    Lognormal    0.079  0.001
## 3        Gamma    0.064  0.012
## 4  Exponential    0.436  0.000
## plot sugar by treatment as a boxplot with violins and jitter
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=sugar_avg, color=Treatment, fill=Treatment))+
  geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  ylab("Sugar Concentration (%)")+
  xlab("Treatment")+
  theme_bw()+
  #facet_wrap(~Treatment)+
  labs(
    title=str_wrap("Tomato Sugar Concentration by Treatment", 40)
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## plot sugar means per plant by treatment as a boxplot with violins and jitter
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=sugar_avg_mean, color=Treatment, fill=Treatment))+
  geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
  scale_fill_manual(values=four_colors)+
  scale_color_manual(values=four_colors)+
  ylab("Sugar Concentration (%)")+
  xlab("Treatment")+
  theme_bw()+
  #facet_wrap(~Treatment)+
  labs(
    title=str_wrap("Mean Tomato Sugar Concentration per Plant across Treatments", 40)
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## scatter plot of sugar by mass grouped and facet wrapped by treatment
ggplot(data=Fl_data_no_BER, aes(x=mass, y=sugar_avg, 
                                fill=Treatment, shape=Treatment,
                                color=Treatment))+
    geom_point()+
    scale_color_manual(values=four_colors)+
    scale_fill_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  ylab("Sugar Concentration (%)")+
  xlab("Tomato Mass (g)")+
  theme_minimal()+
  facet_wrap(~Treatment)+
  labs(
    title="Tomato Sugar Concentration by Mass"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

3.2.3 Sugar analysis

Sugar appears to have either a Normal or Gamma distribution, so we’ll test both.

# Individual Gamma GLM
mod_sug_g_s <- summary(mod_sug_g <- glm(sugar_avg ~ Treatment, data = Fl_data_no_BER,
                                        family = Gamma(link = "log")))
print(mod_sug_g_s)
## 
## Call:
## glm(formula = sugar_avg ~ Treatment, family = Gamma(link = "log"), 
##     data = Fl_data_no_BER)
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.86130    0.01719 108.286  < 2e-16 ***
## TreatmentTransplantation  0.03451    0.02393   1.442 0.149722    
## TreatmentGermination      0.04273    0.02624   1.628 0.103977    
## TreatmentGerm+Trans       0.09020    0.02438   3.699 0.000235 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.04874986)
## 
##     Null deviance: 34.057  on 627  degrees of freedom
## Residual deviance: 33.375  on 624  degrees of freedom
## AIC: 2311.3
## 
## Number of Fisher Scoring iterations: 4
# Grouped Gamma glm
mod_sugG_g_s <- summary(mod_sug_g <- glm(sugar_avg_mean ~ Treatment, data = Fl_means_byplant,
                                        family = Gamma(link = "log")))
print(mod_sugG_g_s)
## 
## Call:
## glm(formula = sugar_avg_mean ~ Treatment, family = Gamma(link = "log"), 
##     data = Fl_means_byplant)
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.85730    0.03448  53.872   <2e-16 ***
## TreatmentTransplantation  0.03351    0.04876   0.687   0.4955    
## TreatmentGermination      0.05741    0.04876   1.178   0.2453    
## TreatmentGerm+Trans       0.09988    0.04876   2.049   0.0465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.0142631)
## 
##     Null deviance: 0.65864  on 47  degrees of freedom
## Residual deviance: 0.59499  on 44  degrees of freedom
## AIC: 117.86
## 
## Number of Fisher Scoring iterations: 4

3.3 Harvest

3.3.1 Harvest exploratory

## scatter plot of ripeness by 'days from harvest to analysis' grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=DaysFromHarvestToAnalysis, y=ripeness, 
                                fill=Treatment, shape=Treatment,
                                color=Treatment))+
    geom_jitter()+
    scale_color_manual(values=four_colors)+
    scale_fill_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  ylab("Ripeness")+
  xlab("Days from Harvest to Analysis")+
  theme_bw()+
  #facet_wrap(~Treatment)+
  labs(
    title="Tomato Ripeness by Days from Harvest to Analysis"
  )+
  theme(
    legend.position="right",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

## scatter plot of sugar by days from harvest to analysis grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=DaysFromHarvestToAnalysis, y=sugar_avg, 
                                fill=Treatment, shape=Treatment,
                                color=Treatment))+
  geom_jitter()+
  scale_fill_manual(values=four_colors)+
  scale_shape_manual(values=four_shapes)+
  scale_color_manual(values=four_colors)+
  ylab("Sugar Concentration (%)")+
  xlab("Days from Harvest to Analysis")+
  theme_bw()+
  #facet_wrap(~Treatment)+
  labs(
    title="Tomato Sugar Concentration by Days from Harvest to Analysis"
  )+
  theme(
    legend.position="right",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

3.3.2 Harvest analysis

do tests n shit

## probability of BER for each mass bin
ggplot(data=Fl_data_mb, aes(x=mass_bin, y=pBER*100, fill=mass_bin))+
  geom_col()+
  ylim(0, 50)+
  scale_fill_manual(values=ten_col)+
  ylab("BER Probability (%)")+
  xlab("Mass Bin (g)")+
  theme_bw()+
  labs(
    title="Tomato BER Probability by Mass Bin"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_col()`).

## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=pBER*100, fill=Treatment))+
  geom_col()+
  scale_fill_manual(values=four_colors)+
  ylab("Percent of Fruit w/ BER")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Tomato Fruit Blossom End Rot Occurrence by Treatment"
  )+
  theme(
    legend.position="none",
    text = element_text(size=30, family="mont", face="bold"),
    axis.title = element_text(size=30, family = "mont", face= "bold"),
    title = element_text(size=34, family="open", face="bold")
  )

## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=fruit_sum, fill=Treatment))+
  geom_col()+
  scale_fill_manual(values=four_colors)+
  ylab("Fruit Count")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Total Fruit Count by Inoculation Treatment"
  )+
  theme(
    legend.position="none",
    text = element_text(size=30, family="mont", face="bold"),
    axis.title = element_text(size=30, family = "mont", face= "bold"),
    title = element_text(size=34, family="open", face="bold")
  )

## Plot blossom end rot by treatment as a stacked bar chart
ggplot(data=Fl_data, aes(x=Treatment, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
  geom_bar(position="stack", stat="identity")+
  scale_fill_manual(values=two_colors)+
  ylab("Fruit Count")+
  xlab("Treatment")+
  theme_bw()+
  guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
  labs(
    title="Tomato Fruit and BER Count by Treatment"
  )+
  theme(
    legend.position="right",
      legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## Plot blossom end rot by mass bin as a stacked bar chart
ggplot(data=Fl_data, aes(x=mass_bin, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
  geom_bar(position="stack", stat="identity")+
  scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
  guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
  ylab("Fruit Count")+
  xlab("Mass Bin")+
  theme_bw()+
  labs(
    title="Tomato Fruit and BER Count by Mass Bin"
  )+
  theme(
    legend.position="inside",
    legend.position.inside=c(.8,.7),
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
    axis.text.x = element_text(size=20, family = "mont", angle=45, 
                               lineheight=1, hjust=1, vjust=1.3),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## Plot blossom end rot by plant as a stacked bar chart
ggplot(data=Fl_data, aes(x=plant, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
  geom_bar(position="stack", stat="identity")+
  scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
  guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
  ylab("Fruit Count")+
  xlab("Plant")+
  facet_wrap(~Treatment)+
  theme_bw()+
  labs(
    title="Tomato Fruit and BER Count by Plant"
  )+
  theme(
    legend.position="right",
    legend.background = element_rect(color=four_colors[3], fill=NA, 
                                     linewidth=.5, linetype = 2),
    text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
    axis.text = element_text(size=16, family = "mont", 
                               lineheight=1),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
    )

## fruit mass sum by treatment
ggplot(data=Fl_data_summary, aes(x=Treatment, y=mass_sum/1000,
                                 fill=Treatment))+
  geom_col()+
  scale_fill_manual(values=four_colors)+
  ylab("Total Fruit Mass (kg)")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Total Tomato Fruit Harvest across Inoculation Treatments"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## fruit mass sum by plant
ggplot(data=Fl_sum_byplant, aes(x=plant, y=mass_sum/1000,
                                 fill=Treatment))+
  geom_col()+
  facet_wrap(~Treatment)+
  scale_fill_manual(values=four_colors)+
  scale_x_discrete(name="Plant", breaks=c("3","6","9"))+
  ylab("Total Fruit Mass (kg)")+
  xlab("Plant")+
  theme_bw()+
  labs(
    title="Total Tomato Fruit Harvest across Individual Plants"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## plot of fruit mass by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=mass_sum/1000, color=Treatment, fill=Treatment))+
  geom_violin(alpha=0.5)+
  geom_boxplot(width=0.2, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
    scale_fill_manual(values=four_colors)+
    scale_color_manual(values=four_colors)+
  ylab("Mass Sum (kg)")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Total Fruit Mass per Plant across Inoculation Treatments"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

## plot of fruit count by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=fruit_sum, color=Treatment, fill=Treatment))+
  geom_violin(alpha=0.5)+
  geom_boxplot(width=0.2, alpha=0.5, outliers = FALSE)+
  geom_jitter(width=0.1)+
    scale_fill_manual(values=four_colors)+
    scale_color_manual(values=four_colors)+
  ylab("Fruit Count")+
  xlab("Treatment")+
  theme_bw()+
  labs(
    title="Tomato Plant Harvested Fruit Count Across Inoculation Treatments"
  )+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )

Fl_avg_sums <- Fl_sum_byplant %>%
  group_by(Treatment) %>%
  summarize_at(vars(mass_sum), list(mean=mean, sd=sd))%>%
  mutate(mean=mean/1000,
         sd=sd/1000)

## mass per treatment columns with error bars
ggplot(data=Fl_avg_sums, aes(x=Treatment, y=mean, fill=Treatment))+
  geom_col()+
  geom_errorbar( aes(ymin=mean-sd, ymax=mean+sd), width=0.4)+
  scale_fill_manual(values=four_colors)+
  ylab("Total Mass (kg)")+
  xlab("Treatment")+
    labs(
    title="Total Tomato Fruit Mass (kg) Across Inoculation Treatments"
  )+
  theme_bw()+
  theme(
    legend.position="none",
    text = element_text(size=24, family="mont", face="bold"),
    axis.title = element_text(size=24, family = "mont", face= "bold"),
    title = element_text(size=30, family="open", face="bold")
  )